{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Variational GP Regression with Pyro"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import pyro\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f61e1f8bc18>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANsAAACXCAYAAACcGS5XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADthJREFUeJztnc1vW8Uaxp9JTCMCunaCVPWW0l67VbfkY3S3CCfcBbtWoUTsEDRRs82lTRUpFZEqaCVLrCqB6D6F3GTHhtb/ACcN64jaagJRVIn0pEiNWtKcu/CxYzvnM8fn+/lJEfaMx3499PW85513niM0TQMhxH+6wjaAkLRAZyMkIOhshAQEnY2QgMj49cZSyiFFUR6a9I0BUAEMKYpy26yNkCThy8ompRwF8KNJ3xAAKIpyH4AqpRwyavPDLkLCxJeVTVGU+1LKikn3xwB+1h9XAIwCeMugzXBVnJmZ4V4FiSxff/21MOvzLYy0IAdgu+n5WyZtpnz55ZemfU+ePMHx48e92OcLtMsdcbTrxo0blmOZICHEAQsLXTh//hjOnDmF8+ePYWHBveuEsbKpAPr1xzkAf+qPjdoMefLkifmbq6pH8/yBdrkjSnYtL/diZqYPu7u1CHF9HbhypRs7O89w4cJzx+8TmLNJKXOKoqgA7gGQenMBwP36SwzaDGlexvf29rC1tYUXL16gXnr27NmzDlreOaJq119//YWenh6cOHECmUwYv7/GRCWMLJWONRytzu5uF0qlfkxOvun4fXyZWT2NL6WUY4qiLOrNDwAMK4ryUNYYBaDWtweM2pywtbWFN954A2+//TaEEPj777/x2muvdfw7eSXKdmUyGaiqiq2tLZw6dSpskyLHxoa7djP8ykYuAlhsaxtuevydwZhDbU548eJFw9HI0RBCIJfLYXt72/7FKeSdd2qho1G7G2KfINE0LdWOpqoqlpaWPL+PEAJpPgFST4C8/vqxQwmQ+fk99Pa2zk1vr4b5+T1XnxF7Z3OK1WQ6ZXZ2Ft9//z1KpRI+/PBDLC0t4ZNPPrEcU6lUUCqVXPc5JZfLYXFx0bS/UjHb7iR1Fha6MDWVwfq6gKYJrK8LTE1lGv9Gxsf3cefOHk6f1iCEhtOnNdy5s4fx8X1Xn5MKZ7ObTKcMDw/j888/R7FYRD6fx8WLFzE2NmY5plAoYHp62nWfG/r6+gzbK5UK7t696/n9k87cXAbPn7dGR8+fC8zNHVxljY/vY23tJR4//h1ray9dOxoQTuo/cKwmc3z8peP3uXjxomHb0tISFhcXkc/ncenSJfzyyy8AgGKxiGq1inK5jPfeew/ffPMNpqenUS6XcfPmTTx48ADlchnFYhGlUqmlT1VVLC4uYmdnB9lsFmNjY8jlco3PVVUVd+/excDAAB4+rOWTKpUKyuVy47N//fVXVKtVrK6uIpvNtvQVCgV3k5hgOpUAsSMVK5vfk1ksFtHX14ebN2/i3XffRX9/P/r7+7G8vIyRkRGoqor3338fuVyu8RxA4/HIyMihvnK5jOHhYeTzeeTz+RZHA2oh7YULFzAyMoJ8Pg+gtlI2f3Z9BR4cHDzURw4wS3S4TYDYkQpnC2Iy685QKpWQzWYxMDAAoHVz1izcM+rL5/OoVqvIZrMYGRlxZIPZZ9evDc3sSjudSoDYkYowcn5+D1NTraHkUSdTVVWUy+VGeDY4OIiVlRWsrq5CVdWGk+zs7KBSqeCnn35CtVpthHSVSqUxFkDjcXtfNpvF4uJiwwnbHW56ehrLy8sYGBhojGn/7KdPn0JVVVSr1UN91WoVg4ODHmY1OdSuv/YwN5fBxkbtR3h+3n0CxA4Rt3TvzMyM1lyI/Ntvv+HcuXON52abxwsLXb5PphVuN7VLpRI+++wz5HI5zM7O4osvvjgUSnbarva5DJO4FiJHreo/FMbH910lQ8JmYGAAKysrAID+/n5fHI0ES2qcLW40h41Or9lItIl9giTtlQ+dIu2VOEEQe2fr6emBqqp0OA9omgZVVdHT0xO2Kb7RiQoir8Q+jDxx4gS2trawvb0NTdOwv7+Prq7o/YZE2a7u7u7GEZskUq8gqmej19eBqakMgGCTZLF3tkwm03IsJI5ZrDCJql2dpFMVRF6J3k8tIR0mqHIsO+hsJPEEVY5lh1+6kWNSylEp5VWDviEppSalfKT/fau339L/O+GHTSS9BFWOZUfHnc2B4Gq/oihCUZSzAD4CcEtvn5BSPkJNN5KQjtGp82he8SNBYiTC2tAU0Z2wjmySQ7jcpFdCSEeJQgWRH2GkI8FVXdznh6amglnoSUgSCDP1/0HzKtd0g40PpJSjbStgC9SN7BxJsmt5uRe3b2exudmNkydf4erVHVe6jn7ZVccPZzMTYW2ncS2nJ0W29TDyT9S0I02x2xeK6r4R7XKHG7sWFrpw/frBftoff2Rw/Xo/stl/dPza7Kjz5UcYeQ8HztIQXJVSNsrWpZTtzqTgQJj1rP6cEMc40REJm447W5Poarvg6oO2l1baxlzSxV0fuRFpJQSIzsa1FX6JtBqJsDaLtFYATNqNIcQpnRJS9RNWkJBEEJWNayvobCQRRGXj2oroXD0S4pEobFxbwZWNkICgsxESEHQ2QgKCzkZiQxR0RLzABAmJBVHREfFCvH4aSGqJQzmWHXQ2EinqoeKZM6daQsU4lGPZEZ+fBZJ4rELFOJRj2cGVjUQGq1AxDuVYdtDZSGSwChXjUI5lB8NIEhnsQsWol2PZwZWNRIYkhIpW0NlIZEhCqGhF4CKtev8hQVa7MSQdjI/vY23tJR4//h1ray8T42hAOCKtQJsgq8MxhMQaP1a2j1FT2AIORFrbuawoytkmuTonYwiJNX5kI52ItBZ0QaAhXS/SkbBrHepGdg7a5Y6o6Uba0i7I6nY8dSM7S5B2LSx0YW4ug42NWkp/ft48AZK0+XIURgoh7rl4T0uRVinlhC5ZBxwIsjoVdiUxpl6Otb4uoGkC6+sCU1OZ2B2VOSpOv+UPQohBIURRCPEvm9faibQaCbIajiHxw+rMWRIq973g1Nl+1jRtFUAfgO+EEP8VQhSNXmgn0mokyGoxhsQIu5UrCZX7XnD6k/JQCLEC4J6maf+pNwohLmqattT+YgcirUb9FGmNOXb3rk5C5b4XnK5s1zRN+7jZsYQQgwD+7Y9ZJI7YrVxJL8eyw5GzaZr2P4O2VU3TZjpvEokrdveuTno5lh3pSAORjmKWBHGyctXLsXZ3XyauHMuOdKSBSMewF97Zc7yPljbobMQVdkmQuJ858xOGkcQVaU/fe4HORlxhlwQh5tDZiCvSnr73Ap2NuCLt6XsvMEFCXMMkyNHgykZIQNDZCAkIOhs5RNxvzRRVeM1GWkjCrZmiCn+ySAtpP+DpJ2HpRk7of7ea2g5pSRL/SPKtmaJK4LqR+mns+/ph0UKT4E+LliTxD6sT1awQ8Y8wdCMLTW0VHGiPtGtJEp9I+q2ZokrgupFt8gdDqIn9AIe1JE2hbqQ3NjZOmbQDxeIWvvqqF7dvZ7G52Y2TJ1/h6tUdFIvPYTHtHSdK89VM7HQjgUa42RD7adeStFrhqBvpDSstkOPHj2NyEpicfAXgld7zpv4XLFGZr3Z81Y10iVMNyFFFUa4BplqSxCcYKoaDH85mpxsJKeVE00o2CmMtSeIBq41pFhOHQ8edzU43Um+/JaV8JKV82jSmRUuy03alCSfKw0m+NVNU8eWazUo3Ur8W63MyhhwNO+kCEg6sIIkpVmEiN6ajCZ0thtiFidyYjiZ0thhiV7/IbGM0obNFGLNQ0S5MZLYxmrCUO6JYHXVxcoMKShdED65sEYX1i8mDzhYiR80oMkyMJwwjQ8LuRLRdqMgwMX5wZQsJZhTTB50tJJhRTB90Np8xkx9wsvGc5nuZJRE6mw1eZN2sKj0YJqYPJkgs8CrrZnVdtrb2ErxxYLrgymaBE1k3LwXBDBPTBZ3NAjtnYUEwcYMvYaR+CFSFiXiPUb/dmDCw2+uyOzc2P7/XEoYCvC5LM2HoRh7qtxsTFnZJDKbviRvC0I006rcb4ytm1112zuImfU/5AeKHs1nqRpr0242xxWw/y8k4q+suqyQG0/fEDbFM/beLtC4v92Jmpg+7uwcp+itXurGz8wwXLjy3fK/Z2X8aXnfNzgoUi9aqpMUiHAuaJlF01E+SaJcfzmanG2nW70RrEsBhkcxS6VjD0ers7nahVOrH5OSbWFjoMt3P2tzsNvyMzc1uR2KcbgRNkyY66jdJsysM3UijfsMxTrFKVDA9T6JC4LqRRv0WYxxh5TCsridRIXDdSIv+I+tGWu1nffqp8VdsTs+zbIoEQSIqSKxS9KyuJ1EhEc4GmO9nMUwkUSExzmYGqzhIVIjlPptbqNdBokAsne3GjRthm0CIa4SmafavIoR4JvHXbIREBTobIQFBZyMkIGKZICHET6SUQ2Ylg15UBmK3skkpx6SUo1LKq0777cYEZNeE/nerqe1WvS9Euw7ZEPZ86af3Nf2+64+klN+a2eqDXaMAfjTp86QyECtni6rkggO7RgHc1+s/C/pzAJiQUj5C7XR6x3H43VtsiMJ8AehXFEUoinIWwEcA6j9Qvs5Xk01m7+9JZSBWzoboSi7YfUahqa2Cg+NElxVFOav/D/YDJ9+93YbQ56ttPqSiKPV//H7Plx2eVAbids0WiuSCV7vaTjQMoXZ+DzhY5fxSFHPy3dttCH2+6uh2/dDU5Pd8+UrcVrZYo4dLjfN7iqLc1n+l32oKLQMlCjZY8IGiKA0dggjYaqQyYKdM0CBuznYUyQXHk+GjXXVGFUW5BjQSJmN6+584CC0Ds8vEhijNV+NaLqD5MqRTKgNxc7bAJRc6ZBeklBNNqeJRAEqTLWf150HbZWRDVOar3ZmCmK96Gl82OTbQIZWB2NVG6mnfCoBC/VpISrlSPwlu0n+oLUi7mtLJ26j9on+kKMp9fcy2PsaXaxCH89ViQ9jzpT8uALimKMpk2xhf58tPYudshMSVuIWRhMQWOhshAUFnIyQg6GyEBASdjZCAoLMREhB0thQihBgTQjwVQuSEED8KIXw7SkMO4D5bShFCjKG2wV7RNC2sKvpUQWdLMUKIFQAjmqZF82ZoCYNhZEoRQowCuIyDg5nEZ+hsKUQIMQHgmqZpDwEUhBB0uABgGElIQHBlIyQg6GyEBASdjZCAoLMREhB0NkICgs5GSEDQ2QgJCDobIQHxf+ow2q5vtkcjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 216x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "train_x = torch.linspace(0., 1., 21)\n",
    "train_y = torch.pow(train_x, 2).mul_(3.7)\n",
    "train_y = train_y.div_(train_y.max())\n",
    "train_y += torch.randn_like(train_y).mul_(0.02)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(3, 2))\n",
    "ax.plot(train_x.numpy(), train_y.numpy(), 'bo')\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "ax.legend(['Training data'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PVGPRegressionModel(gpytorch.models.PyroVariationalGP):\n",
    "    def __init__(self, train_x, train_y, likelihood, name_prefix=\"mixture_gp\"):\n",
    "        # Define all the variational stuff\n",
    "        variational_distribution = gpytorch.variational.CholeskyVariationalDistribution(\n",
    "            num_inducing_points=train_y.numel(),\n",
    "        )\n",
    "        variational_strategy = gpytorch.variational.VariationalStrategy(\n",
    "            self, train_x, variational_distribution\n",
    "        )\n",
    "        \n",
    "        # Standard initializtation\n",
    "        super(PVGPRegressionModel, self).__init__(variational_strategy, likelihood, num_data=train_y.numel())\n",
    "        self.likelihood = likelihood\n",
    "        \n",
    "        # Mean, covar\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(\n",
    "            gpytorch.kernels.MaternKernel(nu=1.5)\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean = self.mean_module(x)  # Returns an n_data vec\n",
    "        covar = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean, covar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = PVGPRegressionModel(train_x, train_y, gpytorch.likelihoods.GaussianLikelihood())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 50/400 - Loss: 24.104   lengthscale: 0.643   noise: 0.474\n",
      "Iter 100/400 - Loss: 19.331   lengthscale: 0.669   noise: 0.310\n",
      "Iter 150/400 - Loss: 14.304   lengthscale: 0.709   noise: 0.193\n",
      "Iter 200/400 - Loss: 6.317   lengthscale: 0.769   noise: 0.111\n",
      "Iter 250/400 - Loss: -0.906   lengthscale: 0.815   noise: 0.056\n",
      "Iter 300/400 - Loss: -9.677   lengthscale: 0.826   noise: 0.026\n",
      "Iter 350/400 - Loss: -15.517   lengthscale: 0.855   noise: 0.013\n",
      "Iter 400/400 - Loss: -19.656   lengthscale: 0.986   noise: 0.007\n",
      "CPU times: user 21.1 s, sys: 28.5 s, total: 49.6 s\n",
      "Wall time: 7.11 s\n"
     ]
    }
   ],
   "source": [
    "from pyro import optim\n",
    "optimizer = optim.Adam({\"lr\": 0.01})\n",
    "\n",
    "def train(num_iter=400):\n",
    "    elbo = pyro.infer.Trace_ELBO(num_particles=256, vectorize_particles=True)\n",
    "    svi = pyro.infer.SVI(model.model, model.guide, optimizer, elbo)\n",
    "    model.train()\n",
    "\n",
    "    for i in range(num_iter):\n",
    "        model.zero_grad()\n",
    "        loss = svi.step(train_x, train_y)\n",
    "        if not (i + 1) % 50:\n",
    "            print('Iter %d/%d - Loss: %.3f   lengthscale: %.3f   noise: %.3f' % (\n",
    "                i + 1, num_iter, loss,\n",
    "                model.covar_module.base_kernel.lengthscale.item(),\n",
    "                model.likelihood.noise.item()\n",
    "            ))\n",
    "        \n",
    "%time train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f61de116320>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAADOCAYAAAAQYJa6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXtwW9d95z8XTwIgCRB8vySRetqyLYmCnDh+JJWUJnabtlHsxHG3q8048bbZbmdnPM1jnY0n2UlbZ9bTTKeZbdO4rdpu6zSKnKRJ/BST2PErgkS9LFmWRIoixTdBEG/g4t6zfwCgSIoviXhcUuczwxGIe3nuTxfAF+f8zvl9jyKEQCKRSPKFqdQBSCSS1YUUFYlEklekqEgkkrwiRUUikeQVS6ku7PP5Ovx+/7F5jj2Wfbje7/d/sYhhSSSSZVKSnorP59sLfH+BY6/4/f7vAO3Z3yUSyQqhJKLi9/tfAbrnOdwO5ISkO/u7RCJZIZRs+DMf2R5Kjg7ge/Od+6UvfUkuspFISsRf/MVfKHM9bzhRyeHz+TqAY/PlXXJ87WtfW7StkZER6urq8hVa3jF6fCBjzAdGjw+WHuOTTz457zHDigqwdylJ2pGRkUUbCgaDeQmoUBg9PpAx5gOjxwf5idEwouLz+Tx+vz+YffyY3+//Zvbx3mwOZk6WqvxG/4YwenwgY8wHRo8Plh9jSUTF5/M9mPnH96Df7z+YffowsDM72/OUz+f7IuAFHlpqu+l0mqGhIZLJJNNrmnRdJxQK5fF/kF8KGZ+iKNjtdhoaGrBYDPMdIlnFlORdlhWSg7Oe25n99xWg6kbaHRoawuVy0dzcjKJczSGpqorVal1GxIWlkPEJIQgGgwwNDdHS0lKQa0gk01lVK2qTySQej2eGoKxkgsEghw4dWlYbiqLg8XhIJpN5ikoiWZhVJSpCiCUJyrPPmti0yYbDYWPTJhvPPnt9t+GJJ57gu9/9Lk8//TQPPPAAhw4d4pFHHlnwb7q7u3n66aev6zoej4eDBw/Oe7y7e76lPjNRFAVpcSFZCF0XvN0T4I1Lyx+G33SD7GefNfH5z1uIxTLic/kyfP7zFiDNww/rS2pj586d7Nu3j66uLrq7u9m3b9+if9Pe3s7jjz9+3fFWVc09Euzu7uaZZ57hG9/4xnW3KZFMZzKucqhrgGOXg2zyLl8SbjpR+epXrwpKjlhM4atftfDww6kltTGXiOzbt49Dhw5x8OBB2tra+OQnP8mRI0cA2L17Nz09PXR2drJ7926efvppHn/8cTo7O+cUhWAwyDPPPMP27ds5diyzTKe7u5vOzs6p9o4fP05PTw9dXV243e4Zx9rb5SJkyeIIITh9JcTBrgE0XVBfaQe0Zbe7qoY/S6Gv7/qevx52795NVVUV3/jGN9i2bRterxev18tzzz3Hnj17CAaD7NmzB4/HM/X7XDzxxBN8/OMfZ8+ePbS1tQGZns709nbv3k1bWxs7duy45phEshjRZJp/P3qFf3q7D6fNnBWU/HDTiUpr6/U9f714PB4Ann76adxuN9u3bwdmLiqab0izEPO1l8vVzHctiWQ2F0ejfKvzIif6J2mpKsNpM+e1/ZtOVL7+9TRO58ykpdMp+PrX09fVTjAYpLOzc2oIAnD06FG6uroIBoO0tbXR09PD8ePH6e7u5mc/+9nUuT09PXR3d8/42+k8/vjjPPfccxw+fHjqnNntTUxMEAwG6enpueZYT0/Pjd8gyaolldb5yclB/va1SyhAo7sMUwFmSpWVPCvwpS99SUyv/blw4QIbNmy45rzZ60CefdbEV79qoa8v00P5+teXnqQtBMVYRzPfvVkqq6lupVSUMr7+iTjf819hNJKkodKO2XStmEzGVdwWjf/+m1sXbe/JJ59ceQWFheThh/UlJ2UlkpWMEILXLwb4yakhnDYzzZ6ygl/zphQVieRmIJXW+eGJQY5cmqCh0o7NUpxshxQViWQVMhFL8f/e7qM/mKClqjC5k/mQoiKRrDJ6x2MceOsyaV0UZbgzm5LN/mRNmOY79qDP59vr8/m+UMyYikE+6nkkkrkQIrPU/v++2oPFpFBbbitJHEY0vu6AqWrl4ELiUyq6urq45ZZbOHToEIcOHeKJJ55Y8t9Or+dZqB4oV9dzIzVDkpsPVdP50YlBDh4boKbcRkVZ6QYhRjS+/hSQW73VzVUTbMOwY8cO2tra2LdvH/v27aOtrY3vfve7S/773OK3+eqBcnU9C50jkeQIxVX+/o1e3uyeoNlThr1ICdn5MGJOxQMEpv1efb0NbH5yXqO46+bc1xbXtJ6eHnbv3j1n7Y/X66WtrY22trZr6nkOHz5MZ2cn+/fv58CBA+zcuXNGm11dXQQCATo7O3n00Ufp7Oykra2NyclJ9u3bx+HDhxetI5Ksbvon4hx46zLxlEazx24I24+bbkVtvpiYmODw4cMcPnyY3bt3s2fPnhm1P0899RQ7d+6cWu06Vz1Prv7nL//yL3n00UfZt28fu3fvnlHXkzvnK1/5Cp/97GfZs2fP1PBpKXVEktXLqf5Jvv3LHoTIFAMaQVDAmD2VIBkbScj0WsYXOnm68bWu66iqyumvfHDGOZqmYTbfWH2DqqpzPi+E4L777ptxnqqqVFRUoKoquq7jcrloa2ujpaWFl19+mXQ6PXUs166u6wghpv7O5XIRDAbRNI333nuPtrY2dF1H13VGR0fxeDxMTExMPXa73de0ORe6ri/JJHw+VoJoGT3GfMUnhMDfF+bFc0G8LgsW3UQ0urzFnCdf9fLy/1tPOGDjW00aX/jCJB//eOyG2jKMqEwzvv4e4Ms+3Q4sOJaZvuw5FArNu9w9n8vgu7q6mJyc5NVXX2XPnj1Tz588eZKTJ08SjUb5sz/7M5555pmpXsef/umf8txzz7F9+3Z6e3s5ffo0AL29vXzrW9/iW9/61lTPZseOHYTDYfr6+ohGo/T29vLXf/3X/OhHP8Lr9fKJT3yC2tpaurq66O3tpa+vb6rNHTt2zBmzyWRa9hJxIy+Bz2H0GJcbn64LXjozwi96E6ytc+dlQduxzkr+42+bUJOZtq5csfDlL3txuytvqHylJLU/WePrvwM+lzO+9vl8R3M+tdm9lLuB9lmbi83gRmt/jIas/ckPRo9xufGpms6Pjg/y9qVMQnau+p0b4Rv7NxAcuXb6ec0awXvvzd0DMlztz0LG19nH8wqJRHIzklQ1nvX3885AmBZPGaY8CQpAcHTuL7Qb9RiSiVqJxOBEEmmeeaOXd4cimSX3eRQUAE/t3Lm4G/UYMkxORSKRXEsgmuIf3rhMIJai0X39MzzHOit5/kAdwVErnlqV+/eP0LF7prn1/ftHOPhXV3MqcGMeQzlWlajkXOONMrVmFOQ9WZkMTSb47uu9qGmdhhuwezzWWTlDLIIjNg7+VRPADGHJPf7ZP9YyOW5lTYuyLI+hVSUqdrudYDC4qvb+WS65zcTs9vx5kEoKz6XxGH//Ri8Wk0JNxY3V8Dx/oG5G7wNATZp4/kDdNb2VZt8om91n2VBt5+8/874bjhtWmag0NDQwNDREIBC4ZttTk8m46aNCxjd921PJyuDMQIh//nUflXYL5cuo4ZkvATv9+VBc5a2eCd4biQKQUDU0XSxrZmlViYrFYplza8/VPtUoWR3ouuC1C+P89NQQNeU2HMs0pPbUqnNOFXtqVRKqxtHLk5zon0QXYFJgS0M5H1jjWvZU9aoSFYlkpTLdpa3RbcdqXn7Pda4ErNWRZtv+d/jnt4dIpjM5k011Lt7fVoUAyizL3/dHiopEUmJCcZV/+XUflwPxvLq05fImmdkfCzW7ruD90LtcUlKQhmZPGR9o907t+TMZn7/M43qQoiKRlJCBYIJ/eLOXeEqj6QamjBejY3eIxp0j/OpigOFQkiRQ5bTygXYv66odBZnQkKIikZSIdwZC/OuRfsosprzuEJhD0wVvdAc40Z/psTitZu5s83BrQ0XeF9BNR4qKRFJkdF3w6vkxfvbOMDWu5Sdk5yKaTPPCmREGJ5OYFNi5xsOO1vwUIC6GFBWJpIik0jrPHR/gaG+QhjwlZGczEEzwwpkRYikNl83MR7fW0egungG2FBWJpEiEE2m+/6tL9E3EaS7AthlCCE5eCfH6xQC6gCZ3GR/dWpf3vZIXQ4qKRFIEesdj/P2vhzHZ7MtOyM5Vz3PbfUF+fm6M86OZRWw7Wt3c1VZV0NzJfJREVLJ+KkGgw+/3f3OB4wv6qUgkRkfXBW92B/jxySGsCtRVLC8hO1c9z3P/XEmX8i5xJYnVrLBncy0b6lz5CP+GKPra9cW24Mj+3p1z3DfiFh0SyVKIpTT+zd/PD08MUldho9y+/GHI7Hoex8Yhah9+k7iSpMpp5ZMdTSUVFCiNn8pStuB4Kvtvu9/vP1aUqCSSPDI4meDbv+jmzECI1qqyvM26TNXtKDqeD75L3b6jmOxpou828FBHE1Wu0mwgNp1SDH8W3ILD7/cf8/l83T6fbwL43GKNLcXM+WYxRC4kMsalIYTg1GCMn70bwG5WcDssxGIZA+lEIrHs9t3VKSJxQc3HuihbG0DoChO/2IKppwn1D06hJm+87XhCw2bVlmWQDgZM1Pp8Pg+ZnsyfA3/n8/mO+f3++TYeW3IhntEL9oweH8gYFyOV1vnpqSHe6I7S6K2gzHrtcMflWt7Q5O7P9HAy3ofZlSIdsTP24x3oI1U8+CcDy247bVJxWLRl38NSiMpiW3A8Bvy53+8P+ny+buBB4JpkrkRiJMYiSf711/0MTCby7iEL2W05eoOcMQUxuyA9WMXgoQ4qnSbu/5OBa/xRSkkpRGXOLTimbdExhd/vP5h11pdIDMvZwRD/dqQfRVFo9uR/kVk8pfHy2VEuT8QB8K31cOcH3Zg+PW8HvqQUXVSyORNfdpP24LRE7GFgp9/v/6bP5/tCtpfilVPKEqOi64JfvDfKC2dGqHbZCrLIbHAywYtnRogkNcosJj58Sy1rq515v04+KdUWHdcIxawtOuRwR2JoEqrGoa4BjvdP0lhpx5Ln5fZCCI73h3izO7M6tqHSzkduraNiGU5wxcL4EUokBiMQTfEvb/cxmM2f5Ms+ILdSdjIEjb93HOvazCTp9pZK7mr35m3zsEIjRUUiuQ4ujcf4p7cuo+mCpjzmT3IrZXGHadx/FKsnjp60cIu1hXs2rAwxySFFRSJZArnZlx90DVBZZqHKmd9tap8/UIe1bYjqB05gsuokhyoZ+1EHYZuVD+++kNdrFRopKhLJIqQ1neffGebV8+PUV9qx59mTRAgBW7qpvesiAJGTLYy/dBtoZoJK8fc6Xy5SVCSSBYgk0zx7pJ/zI5G8boqeI5XWeensKO67YggdJjpvJXx0HZC5znxbkhoZKSoSyTwMhxIceLOPybhKcx4Tsjkm4yo/PTVMIKZiEWaGf9hB5PzV1axWu879+5e3ZH6p6EIQTWq486AIUlQkklmkNR1/7wQ/OTWMzazQ4M6/f2zfRJwX3hkhmdapclr5rdvq6RYpnj+QWnDf43yTTOsEIil0YEt9BR11yx/aSVGRSLIIIegei/HDE4OMhJLUlFvnrN9Z7jVOXgnxqwsBBLDW6+Ajt9Zhs5jo2B0qynJ7IQST8TTRlIbDamb3llo61njwumzLLiYEKSoSCZBZe/LCO8Mc75+kssxCS1X+l9truuDn741xZjACQMcaN+9vq8q7reR8qJrOWERFE4J1Xif7Nlazsa4872bYUlQkNzWptM4bF8d5+d1RFDIbbBXiQx5LabxwLshIRMVsUtizuYZN9eV5v85cxFMagZiK1WziA+u97Frroa4i/3sM5ZCiIrkpEULw7lCYH54YYjKuUltuK9j2Ff0TcV4+O0o0627/wG31BdnnZy7GIil0IXhwRxO3NVfmfTg3F1JUJDcdw6EkPzk1xLnhMFUOa0EqiyFTcHikN8iR3kzxfV25ld+6vQGXvfAfO00XDE4maK5y8MiuFmrKiyNiYFzj6w4ytgj4/f6DRQ5PskrRdcFrF8Z54cwwVpMpr3U7s4kk0rx0dpSByQQISJxo48jLmzlfoxV8Viee0hiNpLh3QzUf3VpflA3EpmM44+ssX86KSbs0vpbkg3AizYG3LvPTU0PUltuorbAVTFB6xqL8m/8KA5MJrMLC2A/uZPjFW0E3ExyxcfCvmjjWWVmQa49HUoSTaf7z+1v52B0NRRcUKE1P5VPAy9nHOePrKXPrbC/mCEgLBEl+uDQe41/e7iOuarRUFa53oumC1y8GOHkl0wtZ63Vw6tt3Er08MyGrJk08f6Aur70VTRcMTSZp8pTxyJ3FHe7MxnDG18AumOrR7JXCIrlRtOxw5/nTQ1Q6LDTkOTk6fVMv7/ogDfu6iCpxTArc1e5le0slr/bN7Rs75YqfB+IpjdFoinvaq7n/tuIPd2Zj1ETteNYhbq/P53twobyKdNMvDistxkhS4ydnAlwYj1PrsmLWdKLRVN6udfJVL//xN42oKTOurf24fvM0UUXDplv5za3l1JZnXPTd1Skmx64VM3d1img0uuw4JmJp0rrgd7ZWs6XORDAwtqz28vE6G9H4epzMsCh37i5gXlGRbvrFY6XE2DMW5V9P9RFXFTY0eAoy3Ol8thVVg+r7T1B+Rz8A0bONTB7dzLrv9k6d98BnRmfsKAiZmp4HPjO6LPd7VdMZDqVoqnHxyK4Wape58+F0VqKb/mLG1wfJOOhDRnSOFD1CyYpE0wW/eG+M508P4y7AcGc6oRDUf/LXlK0JoKsmJl7ZSuRka664eIpc3mT23sc3mk8RQjAWUVF1wUdureWeDTUlH+7MxojG190+ny+YTdhWy5yKZCmE4irfPzFGf0TQ4LZjzbNn7OxrNe1/A7MnSjpsZ+TgLtQRNwCe2muHWLmanmg0uqzeSSylMR5NsbmunN/Z1pjX3kk+Marxde64XKMiWRBdF3T1TfLjk4NEYgnW1roLNrsDVxfOmT066lg5w/9+J1rYARTOqkDTBcOhJE6bmT94Xyu3NVUW9P+4XIyaqJVIFmUknOSHxwe5MBqhptyGQ7EW9MPWMxblxTOjpHVBa1UZjfo6XnaYCUZEQawKhBBMxFRiKZ17NnjZs6WuINuA5BspKpIVR64I8MWzo9hMytTK2DxO7lzDyf4Qr10YRwBbGsr5jU01mE1R7iyQf2xC1RiLpGipcvDo3U20VDkKcp1CIEVFsqK4NB7jB8cGGI0kqauwFTR3ApnewusXAxzvz/RA7lznYdfawswo5a43EklhAvbtaMK3tmrFbM2RQ4qKZEUQTaZ55ewob3SPU1FmKVgR4HTSms7LZ0e5OBbDpMDuzTVsaago2PVyRYAb68p5aGczbkd+HfuLhRQViaERQvDOQIjnjg8ST2k0uvNvPj0X8ZTGT08PMxRKYjObeOC2uoIOQRKqxkgkxYc21vCRW+vyvuNhMZGiIjEsfRNxXjozzLnhCF6nFU8ReicAvYEYv3xvnFAiTYXdzG/f0UC1y1aw603GM8nYR3a1sL2lsLNXxUCKisRwDE4mePnsCO8MhCmzFtaiYDrhRJrXLozTPRYDoLbcxm/fXl8w/xMhBCPhFA6bmc9/sG1FJWMXQoqKxDAMh5J0vjvC8Sshyswmmj2FszycjpZd6+LvDZLWBSZhIupfj/+X67lQXRj/k0z+JEl7rZNP+1qoXKH5k7mQoiIpOeORFD9/bxR/bxCrWaGp0o6pAHmT6VXFuXUl1duHePX8OJPxNADVwsM7z2wnOZ5Z+RocMWf2OIa8CUsyrTMSSnL3hmoeuK2+4DNYxWZJoqIoyveEEJ8qdDCSm4uJWIpX3xvnrUsBzIpCQ6W9YEnY3AboucK+cFzj8MUhHKZhAKqcVu7bWM0//Y/bSY7PzJ/k0/8klFAJJzQe2tmMr4BT06VkqT2Vf1cUZQdQBXQLIS4VLiTJaieSTPPq+XFeuzCGSYH6isKJSY7nD9RlBMWsUbmrB/ddFzDZNIRq5u7NlWxrcWM2KfP6nCzX/0QIwVg0jbvcxh/d18baauey2jMySxWVl4UQIUVRPgF8R1GUl4BjQojOG7noYh610877giwoXD2k0pmd/148M4Ka1qmrsBVt6jQ4asXeHKD6/pNYqzM+JtGzjUz8fAsd/35p6jxPrUpw5NqZnuXsaaxqOkOTSRoqrHz2Q+vxOFdP/mQulvqKHlMU5XuAEEL8phDi/wghOhVF2Xe9F1yiRy3ZKuYPX2/7EuMhhODMQIi/PHyBH58YpMJuptFTVjRB0XRB/UfOUP/Im1iro6jjLoaffR9jP+6gwjHze/X+/SNY7fqM55ZTKDgRUxkJp/it2xv4/Y66VS8osPSeyheFED+Y/kR2OHQncOg6r7mgR61kddE3EecnJwfpGY9ltsMo8rTpWCTFK2dHKduWQugw+eZ6gq9vBM08p1jky/9E0wVDoSS15XYevXstzR5HXrYUXQksSVRmC0r2uS6g6wauuZhHLT6fr8Pv97/i8/m+eAPtSwxAIJri5bOjHOubwGE1F3StyVyzOtt/Y5LjfZO81TOBLqCyzMKa2Dpef7cNdBOeutS8YrHcPY0jyTQTMZV7NlTzkVvqsBdhAy8jYdQpZe/ip0iMSDyl8auL4/z83CgmRaGpsqwg08M5Zs/qBEdsPPePHk4pFwgrmdzJ1sYK7l7vxWbR+FCBqoohM8wbDqUos5r43D3r2FhXnG1NjYbhPGpzvZSlNiaNr4vDQjFmKmtVTg1G6boSRdUFNU4LFrNCPLv+o1D87B/WT/N/FZTf0UfV7jOEFQ2H1cTd6ypo9dhRk3HUZOHiSKZ1xqJpttY7+cgWNy5ijIzEZpyz0l/npWJEj9p2n8/XTkZ4vFmRmTfnIo2vi8fsGCPJNGcHw/zqwjjD4SRmk0KDt6KonqmT2TUlJleC6o+ewrkh8yUTfbeBRx+14SjC0GMsnEI3CT5z7zp2tC5cu7MSX+frxYgetQcBfD7fY2R6MhIDoemCS+Mxfn1pglNXJtEFuMssNLmLs6R+Np7aFKnqMbx738HsVNESFgIv3YZtvAbHH14s6LWFyGzg1egp45FdrVSXF67ocCVhSI/aaedcc56kNEzE05w5N8rr3QHCiTQ2s0JdERatLUQkkWbdf/k1E0omqRrvqWH8+Tswpez87p8MFPTaQggGgkk21pXzyJ0tOFaAzWOxMGqiVmIAhBBcCSZ47fwYb18cxulw4HFai2KQtFhc7wyEeb07gKoIzMJE+PXNjL6xDk9tmvv/cKCgG6Dreua+bG9182BHs+G2yCg1UlQk1yCEoHssRue5US6MRrGbTdSVW6koL62YAARjKp3nxhiYTADQVu3kg5uqKf+NBNHo0WVtgbEUNF0wEEzwgfVePnZH44qzeiwGUlQkU2i64NxwmJfPjDAQSuK0mmjO5kqi0Rtfpp4PdF1wvH+Sty8F0XSBw2rigxtrWF/rLFouR9V0BieTfPiWWvZuqSvoVPlKRoqKhFRa59SVEK+8O0IgqlJZZp4SEyMwGk7SeW6M0UjGLn9zfTn3bPAWZWYnR1LVGI6k+N1tjdy93muYe2NEpKjcpGS2z0xxbjhM57kxoimNKoeFlqrSD3FypDWdI71Bjl2eRAAVdjMf2lzDWm9xK3xjKY1ANMWnfS10rJETkoshReUmQQhBMK7SP5Hg3HCYs4MRYmoaBHhdVqoMVOim64IXf5HmfHwUkyuJENBIDb+zy1X0pGg4kSac1PjMXWvY0lhZ1GuvVKSorGIiiTT9wTjnRyKcGQwTjGXyIlazQqXDisdpnF4JZITv4liMX56eJG5KYnJBariSwMtbGRzz0C4KO6szm2BMJaXp/Nd717FuFfuf5BspKquMQDTFOwMhjl4OMhRKogBmk0Klw0KjgfIks+mbiPNmd4CRcAoUUCecBF/bROxsE5CJOV/ua4uh6RlDaqfNzB/d105TiafQVxpSVFYBSVXj/GiUN7sDXByNogCVjtKtcl2I2RXF9+6/xER1H30TmSlip81M30+2EDmxBvSZQ53luq8thWBMJZLUuHOdh4/cWk95mfyIXC/yjq1QdF3QH4xz7PIkRy8HUTUdp81Mo9uOyWBCkmN6RbGlKoLlrvc4aRqECbCZTXSscbOtpZJv/m3TNYICy3NfW4xUWmcknKK+0s4fvH+NHO4sAykqK4yJWIrTV0K80R0gGFMxmxS8LuuKcGR//kAdmiWF94PnKd/Wh2IS6KoJ9d01fPYPdcqyU8T37x+ZYWcAy3NfW4jc3jsAv3NHA3e2Va2Ie2lkpKgYHF0XDIeTXByNcKI/RN9EHAXwOK0raqyvajpi4yWa7+zOGE7rCuHjrUy+sREtUkbZfz87dW6+3NcWI5xIMxFX2dbs5rdur6fKKQsC84EUFQMSS2lcDsQ4Mxjm9JUQ8bQGAsrLzIbMkyyErgvODoV5+1IQzz0aALHz9Uz8YgvpQMbEyFOXuubvluu+thCqpjMSSuF2Wvnc3evYWOdaUffU6JREVBZz08/aHgCs9/v9q95SMuMYlqB7LMaJ/kl6AzEQYDEreBxWqlzGWUOyVIQQ9AbivHExQCA7lV0uHFw6eAfR7pqp8wo1rJkvpvGISkrX2XtLLfduqL7prB6LQdFFZbqbvs/na59twpT1WXnF7/d3+3y+7/t8vr3X4wRndHQ9swhtJJxkcDLBxbEo568EUKxjKEC53UxDpXGTrUthJJzk9YsBrgQzMzqVZRbe31bFxjoXXSLB8wdSBR3WzEVc1RiLpFhf6+L3tjVRX2kv+DVvVkrRU1nMTb89+/Od7PH2okaXR3RdEIilGA2nuBKM0z0Wo38ijqpltoBQUHDaTVTYzbgrV05+ZC6OdVby4qFylFt6KN+a8TKxW0zsWuvh9ubKqWreQg5r5kIXguFQEqvZxKd2trCj1S0LAQtMKURlQTf9WQZOHWTsJ5eFEGK5TcxLKq0TSqhMxtNMxlVGw0mGw0lGwkkCURWBQCErIDYTVU7LNfvdRPVrcworBVXT6fyl4PSVPioezLysIm0ienwtd3Y42d4aW6SFwhFOpAnGVXauqeKB2+qpkGtOioJh73LhzNwhAAAP50lEQVR2mHRsIX9aWNz4OpzU+L+v9VLhHKLaaaHGZaWm3EqF3Tz1Y5+jnkQIQUoTxFWdmKoTVzViKZ1IUiMYTzMWUwnE0kRTOrnvPQGYFbBZFOxmE+UWZeYwRodkAmb7LycSicVvSImZHmNuC8/zY3G6A0lUk8DeCrpqIvZuE8FfbUQLOXnpRJKt7ztZ9BjTumA0ouJxWPjErV7WVVmIhwLEi9dBmhNpfF04FnTTn8bepSRpFzPpNUWSpJUBHE4HE6rO0FgabUQl81lX0ITAaTVTXW6jymklmkgzmUgTSqikNTF1HkKgZx5hsyjYLFa8lWXUm5W8zBwU2lzoRpm+AraqNcodnzpPqHx0KvkKkBzwEDnZSvRsIyJ1Nak8OW4r6v9LCEFCWEloOg9sa+BDm2oMl4iVxteFYTE3fXw+32O5WaF8JGoVwGE1z+m/IYQgrQvCCZVAJInZrGAzm6h22W56V6/MCthGzE3j1H7gMo4NI1wyCYiBw2pic305tzRW8Defv5VInvcfvl6SqsZgWGVTUwWf2NFEo3tl56hWMoZz088+/1R2d0Iv8FAh41EUBatZyayilBMCU2i64JXX09T8/q+w1UYAELpC7Hwd4nIjf/S1ySnRLeYK2NkIIRiNpBACHtji5cPb2276L4NSYzg3/WyvpKroQUmATOL59ECYE/2TOO/NLFZLh+2Ej64jeroFLVoGisBsupqgKNYK2NkkVI2RcIot9RV8fEcj6WhQCooBMGyiVlJYZlcL79l/BX3tFU4PhEhpmdkybcLFxBsbiJ6ZWeA317CmmFPFuXodRYFP+ZrpaPVgMimMRItyeckiSFG5CZlRLVwdxryzmyPiCkpfRkya3GV0rHEzfqKeH5yfKSjFXAE7F7neyS2NFfzetka8LlmvYzSkqNyEvHCwAvstfVRtHMbRPgqAEKBequPT+8w0ZBfirdsdRmGg6MOaucj1TkwKfHpXC9tb5CI2oyJFZZUxe1hz//4Rtn1oksHJBL2BOL3jMSo/2TN1vkibiJxqIXSknXTQScN/OTujvdywJhqNlmzaO65qjIVT3NpYwe9tb8JjID9dybVIUVlFTB/WmCvipBtHOdw9wlvKKLqiT50nUmbiPTXEu+uIXahHj2WmveaqFi4lGVvHJBazKdM7WWTzc4kxkKKyinjhBxU4dnRTs2UQW/3VIYoOeJ1W1lY7Wet1MNRVy6HnW0oyBbwUctXECU1j19oqPnxLHW6H7J2sFKSorHDiKY0Lo1HODUeofOjqsEZPmUlcqiHeXUu8p5Y//t6lqWMtuyOYDJIrmU0wphJOptlcX8EDt9XLRWwrECkqK4QZuZLGOL7fv0iqZpS+iTh6tl5SqCZi5xuInmkifqkGtMwK4mKbIN0I0WSaiZhKo7uMR+5spb2meNuZSvKLFJUCMFeydDkf4GOdlRz8dj2WxgA17xvAsWGY8yYNApkShLVeB5vqywm9U89zL7UadlgzF0lVYyyiUuGw8MiuVm6bZpMgWZlIUckz05OlAMERGwf/qgnguoQlt6Pg5UCcXw6EaPjDU5isV5Otif4qRF89f/x4DIctW9NUH8Vs0GHNbNKazmg4hcVs4mN3NLBrXVXRdx+UFAYpKnnm+QN1M3oKAGrStKSNsFJpnf6JOL2BOJcn4oQTaQCsLZnjyaFKYucaiZ5pQgs5QRE4vjz3FLBR0fVMrY4m4N6NNdy3sZpyu3wbribkq5ln5tvwaq7nj3S6OPxCGarLSsWmUawNE4hpPf8yq4k1VQ6Of7+N8ZMNU1O/OYpZBbxchBCMR1Mk0oKOVg97t9RSXS5Xw65GpKjkGU+tSnAuG4CGBIOTCUbCSUbCKS4PqcSUJOUPXD1H16ESJ7e22VjjdVJXYUNRFGruNXOwy4o+rT2j50pyCCGYjKeJJNNsqivno7fV0+xxlDosSQExqpv+gseNzP37Rzj47XoUdwRbfQh7YxB74yTWmjA/6Jp2ogJoCsnRClKDHhK9NSQu1VDpFuw/cGFGm6WqAl4ukWSaiViaZo+c0bmZMKKb/oLHjUZS1RiNpBiNpBiLJBkt76fpj0/ArM+OAnhdNuoqbNRV2Hn2f20mNVI5Ne2bIzg6t5+u0XMl00moGuORFG6njf90Zwtbm+SMzs2EEd30FzteElRN561f2nnj5w5S1jjla4KUrw2SVK5dA6IoUOWyUltunxKRmnLbjO00f6y5SGnXOtGtpDzJbGIpjYmYit1i4ne3NbJzrZzRuRkxnJv+Eo7PYCHj6+eec/IXT1UyOHAX7poUex7p5477AnOeK4QgmRaEkhrhpEY4oV19nNSIqzqYwLUHcmV1SUARCtXlZrxOC9VOK9VOS8Yxf8Y3s0YqEWe6/Ox+uI//+Jt1qKmrwmK1aex+uI9o1HjGIPOZc6d1QTCmkRY67jIL97RWcHujC4dVIxgYK2qMRjeWNnp8sHKNr/PKfCa9zz5r4n8+KdCrA7i2quhlKr84kWagTKeyIUEyrZFQdZLpzE9C1aZWps6F0BTSQSfpoBN1wkVqyE1quBKX2cZ/+8fu6477rvuT2O2Dc+RJklyVLWORq1IWQhBKpIkk0ljMJj6wuZadazy0VjlKbkdgdGNpo8cHK9P4ejE3/aW67S/IV79qQbWHaHzQP+P5QWBweO6/0ZMWRMRBndfE2jXgdlhwO6xUOiz870/cPsOsaCpY5cb3FDKCrcD1kFA1AjEVITKreD92ewObGyooM5hjvaS0GNFNf87j10tfH5hcduIXa9ETVrSEFT1pRU9YeeiPRrFbTdgtJrqPVfLiM82kQnZE2gQoDNt1Nv/JAFunJUY9Nem5p4pXcA5kKeRW9o6HUtS4bezZXMu2Fje1FdIlXDI3Rc+i5WZy5nHTX+j4ddHaClqkjJGDdzL2kx1MvHIbk69txtzbwi2NFbTXuGj2OHj1wFqSAScibSY3ZZNbATud+/ePYLXrM55bKWtFbgRNF9ntWpPUVZTx8I5avvzRTey9pU4KimRBDOemP9/x6+XrX0/z+c9biMWujvHnEoGlroBdqWtFrhdV0xmLpBDAHc1u7t1QTbOnjNHR0Wu2a5VI5mLFJ2rn4+GHdSDNE18xc6VfwVM3twjMuwK2xI7xxSauagSiKlazwr0banhfW5U0lZbcEKtWVCAjLHt/O843/uM07Q2eOc8p5UZYpUYIQSSpEUqkcdnNfOz2Bra3unHJAj/JMrjp3z03y7AmhxCCcDJNJKEhgPoKO799ewNbGirkQjVJXrjpRQVW97AGMnYDk4k0sVRmx8HWKgd7NnvYUOei2mWT9TiSvCJFZZWS1nSC8TTJtI4CbKwrZ0erm/ZalzSRlhQUKSqrCFXTCcZUUprAYlbY2lDJttZK1lW7cNrkAjVJcZCissJRNZ2JmEo6KyR3NLvZ1uJmXbVT5kgkJUGKygpE1XQCURVNF1jMJra3eLijpZK1XikkktIjRWUFIIQgoeqEEmk0XWCzmOhY4+GO5krWVjtnWCpIJKVGiopBUTWdUDxNIp0pDahyWvlAu5fNDRWs8TqkkEgMixQVg6ALQSSRJpzUUBSwmk1sqS9na1MlrVUOqpxWOfUrWRFIUSkRqbRONJUmltJJJlTK1CRrq518cFMFbTUuGirLpAWjZEUiRaUIaLogmkwTTekIMv4rZRYz66pdrK91UqbF2NrWfHVTMIlkBWNUN/3Hsg/X+/3+LxY1uDyQTOuEEyrJtMCkgElRWFPl4P11GbuF+ko7HsfV4czIyIgUFMmqwYhu+nuBV/x+f7fP5/u+z+fb6/f7b8ioqVhkRCRNKq0jAJfNwtYmN1vqy2lwl1HtssmhjOSmwYhu+u3Zn+9kj7cXNbpFEEKQ0gThxNUl8E6bha1NlWypL6fZ48DrkklVyc2L4dz0Zxk0dZCxl5yXhdz0AQIxFTWVuiGHek0XxFWdRFonrQtMKOgIym1m2rxlrK9x0Vhhw+MwZ0UkiRZLMhq7vuvcLC7rhcboMRo9PljlbvrZYdKxxewkF3P+NkWSWG1DixpLp9I64WSahKqjKBljSbNJobnGQavXQYvHQbXLSnW5vSB1NDeDy3oxMHqMRo8PDOqmPy3ROp3ubG5kqW75ewuZpNX0zBAmlsqsC7FbzGyur6Ct2klthR2vy4q7zFryLSckkpVGQURlEY/Zxdz08fl8j+VmhfKVqBVCEE1qhJJpEGAyKayvcbG1qYK1Xid1FXYpIBJJHij68Mfv9x/zZZjLTX9n9vmnfD7fF8n0aB5a7jUFgoHJJE2eMt7f7qWtxkWTu0wW30kkBcBwbvrZXklVvq7lddr4/Y46bl3XJL1XJZIisOq/qk0mhbVVZVJQJJIisepFRSKRFBcpKhKJJK9IUZFIJHlFiopEIskrUlQkEkleWfFTIk8++WSpQ5BIJNNQhBCljkEikawi5PBHIpHkFSkqEokkr0hRkUgkeUWKiuQafD7fgz6fb6/P5/vCIucteFxibHLWrvMcW9J7YC5W/OzPdJZgqL3g8WJgdNPvxTyEp523F/gwUPT7uIR72EHWhtTv9x8scni5GJb6XmxfxCqkUPHtBf4WWD/HsSW9B+Zj1fRUpt8IIDhbhRc7bpAYc6bf3wHas78Xm0+RebPDVQ9hw7DE1/HLWTFpN+jr3MFV07LuUsSYu/Y8h5f1Hlg1osLiN8IIH5bFYmif9lypTL8X9BCGzIeihDscLHgPsz2AIwB+v/+b1/MNm0eW8l57Kvtve4liXIhF3wMLsZpEZbEbsawblScWNf2e1hXuAPzFCuw68S5+SsFY7HXcBVT7fL6OEuZ8Fnudj5HpoUzMOm9VsJpEZdWwVNPvArGgh3CJeylLZTx377I9F0Ph8/k8ZO7znwN/5/P5DLUNDUv3kZ6T1SQqi92IZd2oPGEI0+9F+B5Xh10zPIRzz2VnBh4DvCXIByx2D8e5misIkum5FJvFYnwM+PNsAvdzgCGEb9prPOd7YKmsJlFZ7MOwrBuVJxaL8RrT72IHOO0bfi4PYfx+/8FpMyqeOZooNIvdw4PTjnvI5leKzKKvc47svSz6hkDZHpxvVk8u9xrP9x5YEquq9if77dnNtGk6n893NOd/O9dxI8WYfRG/T2ac7QUeWgFDjaKzxNc5AOwqVY9vCTF+IXvcW6r3YqFYVaIikUhKz2oa/kgkEgMgRUUikeQVKSoSiSSvSFGRSCR5RYqKRCLJK1JUJBJJXpGiIikqiqI8qCjKhKIoHkVRvq8oivRkWWXIdSqSoqMoyoNkFvd1CyHk4r5VhhQVSUlQFOUosEcIUfQl6pLCIoc/kqKjKMpeMoV0Ty12rmTlIUVFUlQURXkM+KIQ4hjQriiKFJZVhhz+SCSSvCJ7KhKJJK9IUZFIJHlFiopEIskrUlQkEklekaIikUjyihQViUSSV6SoSCSSvCJFRSKR5JX/D1FBhN02j7+SAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "train_data, = ax.plot(train_x.cpu().numpy(), train_y.cpu().numpy(), 'bo')\n",
    "\n",
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    output = model.likelihood(model(train_x))\n",
    "    \n",
    "mean = output.mean\n",
    "lower, upper = output.confidence_region()\n",
    "line, = ax.plot(train_x.cpu().numpy(), mean.detach().cpu().numpy())\n",
    "ax.fill_between(train_x.cpu().numpy(), lower.detach().cpu().numpy(),\n",
    "                upper.detach().cpu().numpy(), color=line.get_color(), alpha=0.5)\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "ax.legend([train_data, line], ['Train data', 'Prediction'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
